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'  In  a  tokamak,  the  electron  distribution  deviates  from  a  Maxwellian.  Ibis  is  because  the  magnetically 
untrapped  electrons  moving  parallel  to  the  applied  electric  field  tend  to  run  away.  Because  of  the 
presence  of  trapped  electrons,  the  distribution  also  departs  from  the  Chapman  -Enskog  solution  of  tire 
weak  electric  field  problem.  Previous  analytic  and  numerical  methods  have  treated  this  distortion  in  the 
limit  of  vanishingly  small  electric  fields,  a  vanishing  small  number  of  trapped  electrons,  or  both.  We 
present  a  numerical  method  which  relaxes  these  limitations,  and  illustrate  the  distribution  functions 
which  result  from  it.  A 
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NUMERICAL  SIMULATION  OF  TOKAMAK  ELECTRON  DYNAMICS 


I.  Introduction  -  Problems  of  Interest  and  History 

Tokamak  operation  depends  critically  upon  the  physics  of  the  electrons. 
The  electrons  transport  energy,  in  a  way  still  called  "anomalous"  after  a 
decade  of  study.  They  determine  radiation  energy  loss,  they  carry  the 
current  required  for  magnetohydrodynamic  stability.  In  the  absence  of 
collisions,  the  dynamics  of  electrons  trapped  in  local  magnetic  mirrors  is 
qualitatively  different  from  the  dynamics  of  untrapped  electrons.  Electric 
fields  and  collisions  then  give  rise  to  complex  transport  processes.  In 
particular,  an  electric  field  parallel  to  the  magnetic  field  leads  to  a 
parallel  current  density.  For  sufficiently  strong  electric  fields,  this 
current  cannot  be  characterized  by  a  quasi-steady  parallel  conductivity, 
a.  This  is  because  the  runaway  electrons,  which  arise  from  the  decrease 
of  the  Rutherford  cross  section  (with  increasing  energy) ,  make  the  conduc¬ 
tivity  depend  explicity  on  time. 

The  theoretical  determination  of  o  began  with  the  seminal  work  of 

1  2 
Chapman  and  Cowling,  and  was  refined  by  Spitzer,  et  al.  Hinton  and 

Oberman  first  pointed  out  the  non-local  relationship  between  current  den¬ 
sity  and  electric  field  in  toroidal  systems  where  the  mean  free  path  is 

comparable  with  or  larger  than  the  geometric  lengths.  The  more  recent 

U 

developments  are  summarized  by  Hazeltine,  et  al.,  but  all  the  works 
discussed  there  neglect  runaway  electrons. 

The  theory  of  runaway  electrons  was  initiated  by  Dreicer^  and  they 

were  first  observed  in  toroidal  geometry  by  W.  Bernstein,  et  al.**  Further 

7 

analytic  studies  were  carried  out  by  Kruskal  and  I.  B.  Bernstein, 

8  9 

Gurevitch,  and  Lebedev.  References  7-9  were  complex  multiple  domain 
asymptotic  analyses,  the  validity  of  which  was  confirmed  by  the  numerical 
work  of  Killeen1^  and  Kulsrud.11  None  of  the  aforementioned  studies 
included  trapped  electrons. 

Qualitatively,  runaway  behavior  in  even  relatively  weak  parallel  elec¬ 
tric  fields  becomes  important  fof  electrons  of  energy  a  few  times  the  thermal 
energy  or  greater.  The  understanding  of  this  domain  of  energies  was 
Manuscript  rabmittad  March  11, 1981. 
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dramatically  altered  when  Coppi  first  described  what  he  termed  a  "slideaway" 

distribution.  The  first  associated  numerical  calculations  were  carried  out 
13 

by  Hui.  The  implications  of  the  results  as  regards  stability  were  discussed 
lU 

by  Papadopoulos  and  the  resultant  quasi linear  diffusion  further  examined  by 

Hui.^  The  application  of  this  model  to  resistivity  in  tokamaks  was  presented 

by  Wins or . ^  Elaborate  numerical  solutions  of  a  model  Fokker-Planck  equation 

for  times  of  the  order  of  an  electron-electron  collision  time  were  given  by 
17 

Pozolli,  et  al. 

This  paper  is  an  extension  of  the  prior  work  on  runaways.  It  is  not 
concerned  with  toroidal  effects  and  crossed  field  transport,  but  concentrates 
on  phenomena  associated  with  the  parallel  electric  field.  The  model  adopted 
is  that  of  a  very  strong  magnetic  field  with  closed  lines  of  force  along 
which  the  magnetic  field  strength  varies,  and  with  an  externally  applied 
electric  field  parallel  to  it.  The  electrons  are  described  by  the  guiding 
center  kinetic  equation  with  collisions  represented  for 'simplicity  by  the 
Lorentz  form  appropriate  to  electrons  colliding  with  infinitely  massive  ions. 
This  equation  is  presented  in  section  II  where  it  is  written  in  terms  of 
energy,  magnetic  moment,  position  on  the  line  of  force  in  question  and  time. 

The  period  of  the  collision-free  motion  of  a  representative  electron  is  much 
less  than  a  representative  collision  time.  This  permits  reduction  of  the 
problem  to  consideration  of  a  system  of  three  bounce-averaged  kinetic  equa¬ 
tions,  one  for  trapped  electrons,  one  for  circulating  electrons  moving  parallel 
to  the  electric  field,  and  one  for  circulating  electrons  moving  antiparallel 
to  the  electric  field.  The  boundary  conditions  linking  the  equations  are  given. 
The  details  of  the  coefficients  for  a  particular  choice  of  magnetic  field 
strength  are  presented  in  section  III.  A  transformation  of  independent  vari¬ 
ables  which  maps  the  problem  into  a  domain  with  rectangular  boundaries  is 
given  in  section  IV.  This  lends  itself  to  fast  accurate  numerical  solution. 

The  results  are  discussed  in  section  V.  The  details  of  the  numerical  method 
and  the  calculation  of  the  quasi-static  solution  are  presented  in  Appendices 
A  and  B. 


II.  Kinetic  Equation  -  Transformation  to  New  Variables  and  Boundary  Conditions 

The  guiding  center  distribution  function  f(r,  v;  t)  for  electrons 

18 

in  a  magnetic  field  is  governed  by  the  kinetic  equation 
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Here  B  =  ViJixVx  and  E  are  the  magnetic  and  electric  field  vectors  and  b  *  B/|B| . 
In  terms  of  the  velocity  v  one  has  u  =  b  •  v  and  w  *  |b  x  vj ,  S  measures  are 
length  along  B,  and  f=  f(t|/,  x»  s,  u,  v,  t)  where  i|>  and  x  are  flux  coordinates 
which  label  a  line  of  force.  For  simplicity  we  neglect  self  collisions 
amongst  the  electrons  and  represent  the  collision  term  by  the  Lorentz  collision 
operator2 
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which  accounts  accurately  for  small  angle  collisions  of  electrons  with 
ions  and  approximately  for  electron-electron  collisions.  For  electrons 
of  speeds  much  greater  than  the  rms  speed,  Eq.  (2)  is  also  a  good  repre¬ 
sentation  of  electron-electron  collisions.  Here  e,  m  and  n  are  the 
electron  charge,  mass  and  number  density,  Ze  is  the  charge  of  the  single 

species  of  ion  assumed  to  be  present,  to  A  is  the  usual  Coulomb 

2 

logarithm,  X  *  u/v  is  the  cosine  of  the  pitch  angle  and  the  collision 
frequency  is 

„  ,  tojf  t  »  A  .  (I*) 


A  transformation  to  energy  and  magnetic  moment  as  the  Independent 
variables  simplifies  Eq.  (l).  To  see  this,  define 
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and  Eq.  (1)  Is  carried  into 
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where  now  u  ■  +  [H-uB] .  Further  analysis  of  this  equation  requires  some 

knowledge  of  the  orders  of  magnitude  of  its  terms  and  of  the  magnetic  geometry. 

In  order  to  model  the  electron  dynamics  in  a  tokamak,  a  magnetic  geometry 

with  closed  lines  of  force  is  chosen  which  has  a  periodic  variation  in  the  field 

strength  along  the  line  of  force.  One  period  of  the  magnetic  field  is  shown 

in  Fig.  1.  When  collisions  are  negligible  and  b  •  E  is  zero,  an  electron 

with  given  energy  H  and  magnetic  moment  p  moves  on  an  orbit  in  the  (u,s)  phase 

plane  as  indicated  schematically  in  Fig.  2.  Note  that  when  H  <  p  B  the 

max 

electron  is  trapped  between  maxima  in  the  effective  potential  #  «  pB  but  when 

H  >  p  B  the  electron  passes  over  the  peak  in  the  effective  potential,  the 
nkAx 

sign  of  u  never  changes,  and  the  electron  is  said  to  be  untrapped  or  circulating. 
One  can  associate  with  each  orbit  a  quantity 
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-L 
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( untrapped  electron) 


x 


( trapped  electron  ) 


(9) 


For  a  trapped  electron  x  is  one-half  the  conventional  period;  for  an  untrapped 

electron  x  is  the  time  required  for  an  electron  to  move  through  one  space  period 

of  B.  Note  that  as  H  -►  pB  the  characteristic  time  x  *v  -  £n|H  -  pB 

max  max ' 

We  are  concerned  with  situations  where,  for  almost  all  values  H  and  p  and 
for  almost  all  electrons, 

E  ■  lT  “I  *  |vt|  «  1  (10) 


This  inequality  implies  that  the  distribution  function  changes  in  a  time,  deter¬ 
mined  by  collisions,  which  is  much  longer  than  a  representative  period 

time.  Let  us  make  order  of  magnitude  estimates  of  the  various  terms  in  Eq.  (8). 


We  expect  smoothly-varying  solutions  which  behave  in  order  of  magnitude 
like 
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We  seek  a  solution  of  the  form 

f-f  +f  +f  +  ...  (15) 

o  i  2 


where  in  order  of  magnitude  | f ^ ( ^ /f  [ ~e.  Then  on  inserting  Eq.  (15)  in 
Eq.  (8)  and  equating  like  powers  of  e. 
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Equation  (16)  requires  that  fQ  be  independent  of  s;  e.g,,  tQ  m  fQ(H,u,t). 
Henceforth  we  shall  suppress  explicit  indication  of  i|>  and  x-  Thus  if  Eq.  (17)  is 
integrated  one  period  around  an  orbit  there  results 
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The  first  integral  on  the  left  in  (18)  is  t.  We  choose  the  paths  of  integration 
as  in  (9)  for  the  remaining  integrals.  The  reason  for  our  choice  now  becomes 
clear.  The  integral 

D  -  $  yu  ds  (19) 


is  now  continuous  across  the  separatrix. 
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The  annihilation  of  the  f  term  in  Eq.  (17)  was  accomplished  by 
integrating  around  a  period  for  untrapped  electrons  (assuming 
periodicity  of  f^),  and  around  a  closed  orbit  for  trapned 
electrons.  For  untrapped  electrons,  the  integral  of  the  electric  field 
term  has  a  sign  which  depends  on  whether  u  •  E  is  positive  or  negative: 
for  trapped  electrons,  the  corresponding  integral  vanishes. 

Then  Eq.  (18)  implies 
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where  f  denotes  the  distribution  function  of  the  trapped  electrons,  and 
f1  that  of  the  untrapped  electrons,  the  ±  being  taken  according  as 
sgn  u  -  ±  sgn  b  •  E.  These  three  parts  of  the  electron  distribution  are 
indicated  in  Figs.  1  and  3.  On  the  separatrix,  or  better  said  on  adjacent 
points  on  the  side  of  the  item  boundary  layer  straddling  the  separatrix 
in  which  t  ■+  <*> 

f+  =  f~  =  f»  (2U) 

In  addition  f  must  be  regular  at  all  the  natural  singular  points  of 
the  differential  equations.  As  one  approaches  the  separatrix,  where  the 
bounce  average  which  is  used  to  derive  Equations  (20)  -  (22)  fails,  one 
can  invoke  conservation  of  flux  across  the  separatrix  to  provide  the 
required  condition.  The  details  are  given  in  Appendix  A. 


r 


.  3  -  Diagram  in  magnetic  moment  versus  energy  space  snowing  the 
location  of  trapped  and  passing  electrons. 


III.  A  Magnetic  Field  Geometry 

The  transport  coefficients  t,  A  and  D  may  be  calculated  when  the 
magnetic  field  is  known  as  a  function  of  space.  For  simplicity,  we 
shall  assume  the  sinusoidal  form 


B(s)  =  Bmin  il  +  6  Cl  -  cos  (^)]j 
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where 
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This  magnetic  field  behavior  is  depicted  in  Fig.  1.  As  Indicated  in  Fig.  2, 
the  parallel  electron  velocity  varies  with  this  magnetic  field 


u(s )  -  [H  -  nB(s)]!5. 


For  trapped  el«r  rons,  u  vanishes  at  the  turning  point,  ±s0>  where 
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Now  let  us  define 
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Then  in  (H,  p.)  coordinates 


TTS 


max 


2L 


• — r —  26p.B  . 

(•  L  r-  min 

J  1  H-M-B  , 

o  min 


tt[-(H-m.B  .  )•$  o 
mv  ^  min  1 


sin2  9]  d9 


A  - 


.L/2 


2 1  j*  b  •  E  ds|,  H  >  n  B 


max 


,  H  S  (A  B 


max 


(29) 


(50) 


10 


and 


rr8 


max 


*>  =  C  I  (H-^tn5^  -T 


26p.B  . 
[1  - 


sin2  9] 


* 


ttb  .  '  m 

min 


H-|iB  . 

1— ffe  9  -  d9-  CD 


Performing  the  9  integrations  in  Equations (29)  and (31) and  employing  the 
superscript  notation  introduced  in  Section  II,  the  transport  coefficients 
may  be  expressed  by 


t  =  akK(k)  H  >  p.  B, 


where 
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and  K(k)  is  the  complete  elliptic  integral  of  the  first  kind  defined  by 
TT 
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where 
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and  E(k)  is  a  complete  elliptic  integral  of  the  second  kind  defined  by 
TT 

E(k)  =  J*  d9(l-k2sin29)^.  (39) 

o 

For  convenience,  the  9  variation  in  the  denominator  of  the  integral  in 
Equation  (31)  was  neglected,  and  therefore  the  expression  in  Equation  (37) 
is  an  upper  bound  for  the  diffusion  coefficient. 


IV.  Reduction  to  the  Unit  Square  -  Second  Transform  and  Comments  on  Computability 

The  equation  for  the  averaged  distribution  function,  f,  is  in  its 
simplest  form  in  Eq.  (18).  In  velocity  space,  it  consists  of  purely 
convective  "flow"  in  the  H-direction  plus  purely  diffusive  motion  in  the 
^-direction.  This  would  be  easy  to  solve  either  numerically  or  computa- 
tionally,  except  for  the  boundary  and  matching  conditions. 

Extensive  numerical  studies  have  been  performed  on  a  close  relative  of  (13), 
the  Fokker-Planck  Equation.  3*5*11,20  pef6rences  5  and  il  have  not  included  the 
magnetic  trapping  effects,  but  have  required  special  treatment  near  H  *  0 
and  H  =  °°.  References  3  and  20  have  included  trapped-electron  effects,  and 
require  special  treatment  on  the  separatrix.  Much  of  this  difficulty  can 
be  removed  by  choosing  a  compact  coordinate  system  (the  unit  square)  in 
which  the  separatrix  is  parallel  to  one  of  the  coordinate  boundaries. 

This  transformation  may  be  chosen  for  computational  convenience  from 
the  rational  functions.  Consider  the  variables 
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Here  a  and  Hq  are  constants  which  will  be  chosen  later  to  scale  the 
coordinate  system  for  problems  of  interest. 

Equations  (21-25)  when  transformed  assume  the  conservative  form 
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where  fQ  represents  fQ,  f  or  f  with  the  appropriate  definitions  of 
T,  A  and  D  given  in  Section  III. 


13 


1 

v 


(44) 


3x  _  1  1^  (1— x)  (1+ax) 

3H  "  H  y  1+a 

o 


3x  Bmin  1-y  (1+ax)2  (45) 

3p  “  "  H  y  1+a 
o 


-  1_  (1_y)  2 

3H  H  V  y> 
o 

and  Che  Jacobian  Is 


(46) 


Bmln  (1-y) 3  ( 1+ax) 2  (47) 

Hq  2  y  1+a 


The  separatrix  Is  now  along  Che  coordinace  line 

B  -  B  . 

_  max  min 

x  *  x  =  — - : — -  . 

s  B  +  aB  . 
max  min 


(48) 


The  cose  of  chis  simplif icacion  is  a  convecCion  Cerm  wich  derivaClves  in  boch 
coordlnaCes. 

The  numerical  analysis  of  Eq.  (43)  is  now  elemenCary.  The  convecCion 

and  diffusion  operaCions  can  be  performed  on  a  Carcesian  (x,y)  grid  wich 

convencional  mechods.  The  funcCions  f  +  and  f  differ  above  Che 

o  o 

+  - 

separaCrix,  and  are  che  same  (i.e.,  f  -  f  «  f*)  below  ic,  and  are 
conCinuous  on  ic.  The  boundary  condiCions  are 

3f  /3y  -  0  (49) 

o 

aC  x  =  0  and  x  •  1,  and 

3f  / 3x  -  0  (50) 

o 

ac  y  *  0  and  y  ■  1. 

And  Che  flux  conservaCion  condiCioning 

D  ,  3x/ _ o 

J  '■3y;  3x 

conCinuous  on  x  *  x  .  If  desired,  che  conservaCion  condiCion  can  be  relaxed, 
8 

and  an  asympCodc  form  can  be  prescribed  aC  y  ■  1. 

This  form  of  Che  model  is  now  easy  co  apply  Co  physical  problems. 

The  nexC  secCions  presenC  some  appllcaCions.  Appendix  A  describes  the' 
derails  of  che  numerical  methods. 
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V. 


Conclusions 


Figure  1*  shows  the  electron  distribution  function,  /  ,  for  three  values 

e 

of  the  applied  electric  field  at  the  same  final  time.  The  curves  are  lines  of 

constant  /  and  the  difference  between  each  contour  is 
e 

A/  =  exp  [  (*-/i  /  (max)  -  /  (min ) /(N  +l)J  where  N  is  the  total  number  of 

“  c  C  C  C 

contours.  For  case  (a),  the  applied  electric  field  is  so  weak,  E/ED  <<  1,  that 

the  departure  of  the  electron  distribution  function  from  a  Maxwellian  is  not 

discernable  in  the  graph,  i.e.,  lines  of  constant  /  are  semicircles.  In  case 

e 

(b),  the  applied  electric  field  is  larger  and  the  departure  from  Maxwellian, 

though  small,  is  row  evident  in  the  plot.  Finally,  in  case  (c),  the  electric 

field  has  generated  a  large  high  energy  tail  on  the  distribution  function  and 

the  level  curves  are  markedly  distorted  from  semicircles. 

Figures  b  and  6  display  the  electron  distribution  function  as  a 

function  of  the  parallel  and  perpendicular  velocities  respectively,  i.e., 

/  (vi)  and  /  (v  ).  For  the  small  electric  field  case,  the  distribution  func- 
e  it  ex 

tions  remain  essentially  Maxwellian  with  equal  parallel  and  perpendicular  temper¬ 
atures  (the  initial  condition).  For  the  case  of  the  large  electric  field, 

/@(vj)  is  highly  asymmetric  and  the  increase  in  the  perpendicular  temperature 
is  attributable  to  pitch-angle  scattering. 

Figure  7  shows  the  variation  of  the  plasma  resistivity,  n,  defined  as 
the  ratio  of  the  current  density  to  the  parallel  electric  field,  as  a  function 
of  time  for  two  values  of  the  applied  electric  field.  Each  curve  exhibits 
two  different  types  of  behavior.  The  first  is  a  very  rapid  response  to  the 
electric  field  which  is  shown  by  the  sharp  decrease  in  the  plasma  resistivity. 
However,  as  the  distribution  distorts  from  the  initial  Maxwellian,  the  pitch- 
angle  diffusion  becomes  important  and  a  second  regime  is  reached  with  a  much 
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Fig.  4  -Velocity  apace  plots  of  the  electron  distribution  function  for 
their  values  of  the  applied  electric  field  after  fifty  colllaion  tinea 
a)  E/Ed  -  .0085,  b)  E/Ed  -  .085,  and  c)  E/Ed  -  .85. 
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Tig.  5  -  The  electron  distribution  function  es  e  function  of  the  parallel  velocity 
The  aolid  line  represents  the  initial  distribution  end  the  dashed  line  represents 
the  distribution  after  fifty  collision  times. 


T 
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Fig.  6  -  The  electron  distribution  function  as  a  function  of  the  perpendicular 
velocity.  The  solid  line  represents  the  initial  distribution  and  the  dashed 
line  represents  the  distribution  after  fifty  collision  times. 


v'v 


Fig.  7  -  The  resistivity  as  a  function  of  time  for  two  values  of 
the  applied  electric  field.  In  the  upper  curve  E/Ed  *  .0085  and 
in  the  lower  curve  E/Ed  ■  .85. 
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more  slowly  varying  resistivity.  The  departure  from  a  constant  resistivity 
is  attributed  to  the  runaway  electrons  and  can  be  seen  to  be  present 
even  for  the  small  electric  field  case,  though  it  sets  in  at  ever  later 
times  as  the  electric  fie±d  decreases. 

Figures  8  and  9  display  the  results  of  varying  the  number  of  trapped 
particles.  Figure  8  shows  the  electron  distribution  function  for  three 
values  of  6  =  (e  =  a/R)  at  the  same  final  time  for  the  same  applied  electric 
field.  The  variation  in  the  number  of  trapped  particles  is  shown  by  the  size 
of  the  angle  subtended  by  the  two  dashed  lines  representing  the  position  of 
the  separatrix.  The  time  evolution  of  the  plasma  resistivity  corresponding 
to  these  cases  are  shown  in  Fig.  9-  The  important  point  to  note  here  is  the 
increase  in  the  resistivity  due  to  the  decrease  in  the  available  current 
carriers,  i.e.  ,  trapped  particles  cannot  contribute  to  the  plasma  current. 

Finally,  Fig.  10,  shows  the  time  evolution  of  the  plasma  resistivity 
for  three  different  initial  electron  distribution  functions.  The  solid  line 
represents  the  resistivity  resulting  from  an  initially  Maxwellian  distribution. 
The  dashed  line  represents  the  resistivity  when  the  initial  distribution  was 
that  calculated  in  Appendix  B  and  corresponds  to  the  counterpart  of  the 
Spitzer-Harm  calculation  for  the  case  of  trapped  particles.  The  ratio  of  the 
two  asymptotic  values  was  found  to  be  approximately  5/3.  When  this  factor 
was  included,  and  it  is  believed  to  be  a  numerical  error  in  the  code,  the 
result  was  the  dotted  line.  The  initial  transient  is  believed  to  be 
associated  with  finite  grid  size  effects.  This  shows  that  the  Maxwellian 
plus  the  first  order  correction  yield  a  much  better  approximation  to  the 
asmyptotic  distribution  than  the  Maxwellian  alone. 
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8  -  Velocity  apace  plota  of  the  electron  di8tribution  function  for 
three  values  of  6 ;  a)  6  *  .100,  b)  6  *  .175  and  c)  6  *  .250 
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In  summary,  these  results  confirm  the  intuitive  notions  that  trapped 
particle  populations  inhibit  conduction,  that  the  Dreicer  field  roughly 
separates  situations  of  approximate  quasi-static  solutions  from  those  of 
well  developed  runaway,  and  that  a  Spitzer-Harm  type  of  solution  is  a  good 
approximation  for  suitably  short  times  in  the  weak  field  case. 
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Appendix  A  -  Details  of  the  Numerical  Methods 


Equation  (28)  with  the  prescribed  boundary  values  and  suitable 
initial  data  constitutes  an  initial-value  problem  for  the  averaged 
electron  distribution  function  F(x,y).  This  section  describes  the 
resulting  implicit  difference  equations,  and  an  algorithm  for  their 
solution. 

Equation  (28)  has  the  form 


6G 

St 


b (AG  > +  h  Cb  w  +  06 ]' 


(Al) 


(.1  (.1  +  ■■  . 

Thus  j  j  (G  ^  +G'  dxdy  changes  only  if  the  boundary  conditions  allow  a 
flux  through  the  boundary.  This  property  will  be  preserved  in  the  numeri¬ 
cal  algorithm.  It  can  most  easily  be  done  if  the  x  and  y  derivatives  in 
Eq.  (Al)  are  replaced  by  centered  differences. 

Let  the  boundaries  of  the  Cartesian  grid  be  x  =  I/N,^  1*0.. N^  and 
y  =  J/Ny,  J  =  O..Ny.  Then  the  change  in  G  within  a  grid  cell  is 


g  (Hi,  j4)  -  Cp  (i+i,  j4)  -  p  (i,  j4)]  +  ^  [q  ini,  J+i) 

-  Q  (Hi,  J)]  (A2) 

in  terms  of  the  fluxes  through  the  cell  boundaries 


p  (i,  j4)  *  i  i  [c(i+i,  J-*i)  g  (i+i,  J+i)  +  c(i-i,  j +i)  g  (i-i,  j-i)] 

.  i_  [B(i+i,  j+i)  g  (i+i,  J+i)  -  B(i-fe,  J+i)  g  (i-i,  j+i) 

6x 

+  B(i-i,  J+i)  g  (i+i,  M)  -  B(i+i,  j+i)  g  (i-§,  J+i)]}  (A3) 

and 

q  (ni,  j)  =  itA(ni,  j-ti)  g  (nl,  j+i)  +  A(ni,  j-i)  g  a+i,  J-i)> 

(Ak) 

These  expressions  apply  on  the  Interior  boundaries  of  all  cells.  On  the 
edges  of  the  unit  square,  the  fluxes  are  zero.  On  the  separatrix,  the  two 
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fluxes,  P+  and  P-,  both  represent  transport  from  the  density  in  the 
adjacent  trapped  electron  cell,  while  they  add  separately  to  their 
respective  untrapped-electron  densities. 

The  separatrix  is  most  conveniently  treated  by  choosing  a  i.i  Eq.  (10) 
so  that  it  lias  on  a  grid  boundary.  In  particular,  for  a  physics  problem 
in  which  both  trapped  and  untrapped  electrons  are  important,  and  there  are 
no  external  sources  of  electrons  at  some  particular  energy,  the  most  natural 
choice  of  the  adjustable  grid  parameters  is 

a  =  ^2  -2  (A5) 

min 

and 

Hq  =  kT,  (A 6) 

where  T  is  a  representative  temperature  of  interest;  e.g.,  the  anticipated 

thermal  electron  temperature  of  the  plasma.  Cases  with  B  <2B  ,  cause 

max  min 

no  problem,  since  a  appears  in  Eqs.  (25-32)  only  in  the  combinations  (1+a) 
and  (1+ax),  Thus,  no  poles  or  zeros  are  introduced  in  the  unit  square 
provided  a  >  -1. 

It  remains  to  specify  the  time  difference  formulation  of  Eqs.  (A2-A^)- 
It  can  be  most  compactly  stated  in  terms  of  Eq.  (Al).  A  convenient  general 
form  is 

[G  (t+1)  -G  (t)]  =  ( 1— S )  ~  [AG  (t)3  +  { 1— S )  [B  +0G  (t)] 


+  B  ^  [AG  (t+1)]  +  s  ~  [B  fCG  (t+1)],  (AT) 

where  G  (t+1)  represents  evaluation  of  G  at  time  t  +  6t.  Here  s  *  1  is  an 

explicit  algorithm,  s  *  0  is  fully  implicit,  and  s  *  0.5  is  an  implicit 

algorithm  which  is  second-order  accurate  in  time. 

Suppose  the  parameters  in  Eqs.  (A5,  A6)  are  used.  Then  G  consists  of 

I  N  Jj  elements,  of  which  N  N  represent  G+,  4  N  N  represent  G~,  and 
x  y  x  y  x  y 

i  N  N  represent  G*.  A  grid  of  this  type  is  shown  in  Fig.  11.  Treating  all 
x  y 

these  elements  as  a  column  vector,  Eq.  (AT)  with  the  space  differences  expanded 
according  to  Eqs.  (A2-AU),  leads  to  a  matrix  equation  of  the  form 
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(Jl  I  -  8  L  )  G  (t+1)  -  (“  I  +  (l-s)  L  )  G  (t).  (A8) 

Here  L  is  the  space-difference  operator  and  I  is  the  identity  matrix. 
Now  the  right-hand  side  can  be  explicitly  calculated,  and  the  left-hand 
side  can  be  solved  by  any  conventional  sparse-matrix  equation  solver. 
This  is  all  that  is  needed  to  advance  the  solution  of  Eq.  (13)  in  time, 
and  the  distribution  function  can  be  calculated  from  it  by  retracing  the 
coordinate  transformations. 
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Appendix  B  -  First  Order  Correction  to  the  Electron  Distribution  Function 

We  desire  a  quasi-steady  solution  of  Eqs.  (20-22)  for  which  the  collision 
term  is  formally  dominant.  Seeking  a  solution  of  the  form  given  in  equation  15 
we  find  the  first  order  equations  to  be 

uher*  i  •  -•  *•  B1 

The  solution  of  Eq.  (IB)  satisfying  reaularity  at  U  =  0  and  vanishing  at 
inifinity  is 


The  next  order  equations  can  be  written 


For  trapped  electrons  the  solution  is  again 


B2 


B3 


'i-'i (H) 


and  can  be  absorbed  into  f  For  circulating 
integrated  to  obtain 


electrons , 


Eq.  ( 2A)  mav  be 
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y  +  c . 


D<  ±  A  < 
D  3  y  "  ±  A  3  H 


BU 


Since  D  vanishes  as  y  -*■  0  the  constant  of  integration  must  be  zero.  A  further 
integration  yields 


-  *  * 


<  f 

3H  J  D(y-) 


+  c 


B5 


In  order  to  make  /  continuous  on  the  separatrix  the  constant  of  integration  is 

absorbed  into  the  integral  form  by  writing  it  as  a  definite  integral,  i.e., 

H/B,. 


MAX 

if? 


b6 


Finally  we  have 


H/B, 


3f 


'+='o(H)+A  -3H 


3/ 

3H 


is/ 


f-7 


MAX 


dy^y' 

D(y') 


i/BMKX 


D(y')  ’ 


B? 


B8 


*  # 

f  *  /(H)  , 


B9 


+  -  * 

and  we  have  taken  /  (H)  a  /  (H)  =  /  (H)  to  be  Maxwellian, 
o  o  o 
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